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FAST  FOURIER  TRANSFORM  IMPLEMENTATION 
FOR  THE  CALCULATION  OF  NETWORK  FREQUENCY 
DOMAIN  TRANSFER  FUNCTIONS  FROM  TIME  DOMAIN  WAVEFORMS 


by 

William  L.  Cans  and  N.  S.  Nahman 


ABSTRACT 

This  report  is  concerned  with  the  software  applications  of  the  fast 
Fourier  transform  algorithm  to  the  relationship  between  time  domain 
waveforms  and  frequency  domain  spectra.    The  first  chapter  is  devoted 
to  a  description  of  the  discrete  Fourier  transform  and  the  fast  Fourier 
transform.    Chapter  2  contains  the  text  and  a  brief  description  of  all 
FORTRAN  n  programs  utilized  in  connection  with  this  work.    All  com- 
putation was  performed  on  the  in-house  time  share  computing  system  in 
the  NBS  facilities,  Boulder,  Colorado.    In  Chapter  3,  problems  en- 
countered using  the  fast  Fourier  transform  algorithm  are  discussed,  an 
example  of  a  time  domain  to  frequency  domain  calculation  is  presented, 
and  future  developmental  considerations  are  mentioned.    In  addition 
Appendix  A  contains  a  detailed  example  aimed  at  disclosing  the  inner 
mechanisms  of  the  fast  Fourier  transform  algorithm. 


Key  Words:         Discrete  Fourier  transform;  Fast  Fourier  transform; 
Frequency  spectra,  discrete;    Network  transfer  function;  Time  domain 
waveform;    Transfer  function. 


1.    DESCRIPTION  OF  FAST  FOURIER  TRANSFORM 
1.  1     Definition  of  DFT. 

For  many  years,  the  Fourier  integral  transform,  the  Fourier 
series,  and  the  Laplace  transform  have  been  available  to  engineers  for 
the  purposes  of  time  domain- frequency  domain  correlations.  More 
recently,  the  sampling  theorem  and  the  discrete  Fourier  Transform 
(DFT)  have  evolved  as  useful  tools  for  handling  sampled-data  band- 
limited  time  waveforms  and  frequency  spectra.     Thus,  while  the  Fourier 
series  and  the  Fourier  and  Laplace  transforms  are  well  suited  to  opera- 
tions involving  continuous,  or  piece-wise  continuous  functions,  the  DFT 
allows  operations  to  be  performed  directly  upon  a  series  of  discrete 
data  samples  representing  a  continuous  time  fvinction. 

Most  recently,  (1965)  a  method  has  been  reported  [  1]  whereby 
the  Fourier  coefficients  of  the  DFT  may  be  machine  calculated  in  a 
fraction  of  the  time  previously  required.     This  algorithm  is  commonly 
referred  to  as  the  Fast  Fourier  Transform,  or  FFT. 

In  order  to  discuss  the  FFT  it  is  first  necessary  to  define  the  DFT. 
Definitions  vary  in  the  literature  ,  however,  this  paper  will  follow  the 
definitions  and  notations  used  in  [  2]  .    The  DFT  is  defined  by 


where  Ar  is  the   rth  coefficient  of  the  DFT  and  X^^  is  the  k  sample 
of  the  time  series  being  transformed.     j  =/^l     and  there  are  N 
samples  in  the  total  time  series.    Thus,  a  time  waveform  consisting  of 
N  sequential  samples  will  yield  a  transform  consisting  of  N  frequency 
domain  coefficients.    For  convenience,  let 


N-1 


Ar 


r  =  0,1,2.  .., 


N-1 


(1.1) 


th 


W 


-27r  j/N 


(1.2) 


2 


then,  in  simplified  form, 

N-1 


'  k=o 


As  in  the  case  of  the  Fourier  integral  pair  there  exists  an  inverse 
transform  of  the  DFT.     This  transform  is  written 

N-1 

X.  =    ^m^^W"''^  i  =  0,l,2...,    N-1  (1.4) 

i       N  ^ 

r=o 

and  is  called  the  inverse  discrete  Fourier  transform  or  IDFT. 

1.  2   Some  Properties  of  the  DFT 

Both  the  Ar's  and  X^'s   may  be  complex  numbers;    the  Ar's,  in 
general,  are  almost  always  complex,  but  for  a  real  time  series  wave  - 

form,  the  X  's  will  be  real. 

i 

In  addition  there  exists  a  convolution  relationship  between  the  DFT 
and  IDFT,     The  IDFT  of  the  product  of  two  DFT's  is  the  periodic  mean 
convolution  of  the  two  time  series  waveforms  of  the  DFT's. 

Likewise,  most  other  properties  of  the  Fourier  integral  trans- 
form may  be  shown  to  apply  to  the  DFT.    Asa  further  example,  the 
DFT  of  the  sum  of  two  time  series  waveforms  is  the  s\im  of  the  DFT's 
of  the  two  waveforms. 

Thus,  with  few  exceptions,  ftie  DFT  lends  itself  to  time- frequency 

operations  on  discrete  sampled-data  functions  with  the  same  power  and 

versatility  as  the  Fourier  integral  transform  brings  to  operations  on 

continuous  functions.    One  exception  to  be  kept  in  mind,  however,  is  that 

data  obtained  by  use  of  the  DFT  is  valid  only  at  the  points  in  time  or 

frequency  where  samples  were  taken.    Interpolation  between  the  Ar's  or 

X  's  is  neither  valid  nor  wise. 
i 
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1.  3    Definition  of  the  Fast  Fourier  Transform 

The  fast  Fourier  transform  is  an  algorithm  which  allows  the 

machine  computation  of  a  set  of  DFT  coefficients  to  be  accomplished 

much  faster  and  more  efficiently.       Previous  straightforward 

2 

methods  require  approximately     N       arithmetic  operations  to 

transform  a  time  series  waveform  consisting  of  N  =  2^   samples.  By 

use  of  the  FFT,  the  same  results  can  be  obtained  with  only  2N  log  N 

operations.    A  1 024  point  waveform,  for  example,  would  require  about 
2 

(1024)    or  1,  048,  576  arithmetic  operations  to  obtain  its  1024  DFT 
coefficients  by  conventional  means.    Using  the  FFT,  these  same  co- 
efficients may  be  obtained  after  only  (2)  (1024X10)  or  20,  480  operations. 
Thus,  the  time  and  effort  required  to  transform  a  1024  point  time  series 
using  the  FFT  is  about  l/50th  of  that  required  using  straight-forward 
techniques.    With  larger  time  series,  the  savings  increase. 

Although  many  variations  of  the  FFT  algorithm  now  exist,  only 
Cooley  and  Ttikey's  decimation  in  time  algorithm  shall  be  described.    [  1] 
Suppose  we  wish  to  calculate  the  DFT  of  a  time  series  waveform  consis- 
ting of  N  samples.    For  simplicity,  let  N  =  2''^  where   n  is  an  integer. 
The  DFT  of  this  waveform  is  defined  as: 

A^   =5]^X^e-2'^j"^/^  r  =  0,   1.  2...,    N-1  (1.5) 

^  k=o 

If  the  original  time  sequence  is  now  split  into  two  sequences  such  that 

the  even-nxrmbered  X^^'s  become  one  sequence  and  the  odd-nximbered 

Xj^'s   the  other  sequence,  then  we  may  define, 

=  X,, 
k  2k 

N 

and  k   =   0,  1,  2,  .  .  .  ,    -J  -1 
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N 

Both  Y     and   Z,    are  now    -—  -point  time  series  and  they  each  have 
k  k  2 

DFT's  as  follows, 


and 


B  = 


:=o 
N 


-4  ujrk/N 


C  Z,  e 

r      4^  k 

k=o 


-4  7rjrk/N 


N 

r  =  0,  1,  2,  .  .  .  ,    —  -1 


(1.7) 

N 

r  =  0,  1,  2,  .  .        -  -1 


Now,  writing  the  original  DFT  in  terms  of  these  two  new  ones: 

r  =  0,  1,  2,  .  .  .  ,  N-1 


or 


k=o 


k  k 


^         Y    e-4  7rjrk/N  ^  ^-2  TTjr/N  ^-4  TTjrk/N 

r      Z-»     k  ^        k  ^ 

k=o 


using  (1.7), 


(1.8) 


(1.9) 


A    =B  +e-^^j^/^C 


r  r 


N 

r  =  0,  1,  2.  .  .  ,    —  -  1 


(1.  10) 


N 

For  values  of   r   greater  than    -—  ,  the  DFT's  of  B    and   C  repeat 

2  r  r 

N  N 
periodically  their  values  for    r  <  -r  .     Therefore,  substituting    r  +  -r 

for   r  in  ( 1.  10)  yields 


5 


=  B    +e    N    '•        Z-'  C 
r  r 


0  ^r<  - 


(1.  11) 


or,  since  e 


-J  77 


1  . 


-  2  ffjr/N  _ 

e  Li 


N 

0^r<  - 


(1.  12) 


Therefore  (1.  2).  (1.  10)  and  (1.  12)  may  be  combined  to  yield        in  terms  of 

B     and   C  : 

r  r 


A    =  B    +  W  C 
r        r  r 


(1. 13) 


A(r+^1    =   B    -  W^C  0 
\         2;  r 


N 


(1.  14) 


Thus,  the  DFT  of  X     may  be  obtained  by  computing  two  smaller  DFT's 
N 

of  length  —  .    It  should  be  noted  that 


(1.  15) 


If  this  decimation  process  is  continued,  B    and  C    would  each  reduce  to 

r  r 

N 

two  —  -point  transforms.     The  resulting  four  transforms  would  reduce 
N 

to  eight  —  -point  transforms,  and  eventually,  the  problem  would  be 

^  N  n 

reduced  to  solving  — -  two-point  transforms.    In  particular,  for  a   2  - 

point  waveforn,  n   such  reductions  will  be  necessary  to  reduce  the 

problem  to  one  of  solving  only  two-point  transforms. 

The  result  of  this  reduction  process  is  a  problem  consisting  only 

of  complex  multiplications  and  additions.    Furthermore,  three-fourths 

of  the  multiplications  may  be  eliminated  because  they  either  involve 

multiplication  by  unity,  or  they  are  red\indant  multiplications.  Thus, 


6 


in  general,    an  N-point  transform  will  consist  of  N  log  N  coraplex 
1 

additions  and  —  N  log^  N  complex  multiplications.  See  Appendix  A 
for  a  detailed  example  of  an  eight-point  transform. 
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2.    OPERATING  SYSTEM  SOFTWARE 

2.  1    FFT  Subroutines 

Four  subroutines  are  used  to  compute  the  DFT  and  IDFT 

14 

of  a   time   waveform   consisting  of  up  to   2     ,    or    16,  384  sequential 
data  points.    These  subroutines,  as  well  as  all  other  programs  in  this 
report,  are  written  in  FORTRAN  II  and  are  implemented  on  the  in-hou 
time  share  computing  system  at  NBS,  Boulder,  Colorado. 

In  order  to  compute  a  "forward"  DFT,  or  to  go  from  the  time 
domain  to  the  frequency  domain,  the  following  call  sequence  of  sub- 
routines is  used: 

CALL  REORDER  (A,  B,  MM) 

CALL  CFFTRC  (A,  B,  MM,   .  5>:^SC,  NX) 

CALL  REALTRAN  (A,  B,  MM,  NX,  INV). 

The  notation  used  in  these  call  statements  relates  to  the  DFT  as 
follows: 


N-1 

S       =   A    +jB     =  SCV  X  w|^^^^^^ 
r  r  r  k  N 

k=o 

where  NX   =   ±1  (indicates  exponent  sign  of  W^^) 

SC    =    real  scaling  factor 

INV    =   ±  1  (indicates  forward  or  reverse  transform) 

^(IvIM)+l  ,         ,  r 

2  =   length  of  transform 


As  an  example,  for  the  case  of  a  1,  024  point  DFT  the  specific  call 
statements  would  be: 
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CALL  REORDER  (A,  B,  9) 
CALL  CFFTRC  (A,  B,  9,  .S^H.,  1) 
CALL  REALTRAN  (A,  B,  9,   1,  1) 
The  se  quence  of  call  statements  for  the  IDFT  to  return  from  the 
frequency  domain  to  the  original  time  domain  waveform  are: 

CALL  REALTRAN  (A,  B,  MM,  NX,  INV) 
CALL  CFFTS  (A,  B,  MM,   l./(N-SC),  NX) 
CALL  REORDER  (A,  B,  MM) 

For  the  previous  example,  the  sequence  of  call  statements  for  the  IDFT 
relating  to  the  1024  point  original  time  series  would  be: 

CALL  REALTRAN  (A,  B,  9,   -1,  -1) 
CALL  CFFTS  (A,  B,  9,   1./1024.  ,  -1) 
CALL  REORDER  (A,  B,  9) 

For  the  1024  point  case,  SC  is  assumed  to  be  unity;  NX  and  INV  are  +1 
for  the  DFT  and  -1  for  the  IDFT. 

A  source  of  possible  notational  confusion  exists  with  the  two 
arrays  A  and  B  since  they  are  used  to  define  both  input  and  output  data. 
For  the  DFT,  or  time  domain  to  frequency  domain  case,  A  and  B  con- 
tain the  N  values  of  time  domain  data  points  as  follows: 

N 

0,   1,  2.  .  .  ,  j-l  (2.  2) 

0,   1,  2.  .  .  ,  j-1  (2.  3) 

Thus  the  first  half  of  the  time  domain  waveform  is  contained  in  A 
and  the  second  half  in  B. 
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Once  the  DFT  subroutines  have  run,  the  contents  of  A  and  B  are 
replaced  with  the  output  data  as  follows: 


S     =  A  +jB 


r  =  0,  1.  2,  --1 


(2.4) 


r  r  r 


where   A     and   B     are  the  real  and  imaginary  parts,   respectively,  of 
r  r 

the  coraplex  frequency  domain  coefficients.    For  a  1024  point  time  do- 
main waveform,  then,  the  first  512  points  would  be  sequentially  stored 
in  A  and  the  last  512  points  in  B.    The  resulting  DFT  solution  would  be 
512  frequency  coefficients  each  consisting  of  a  real  part.  A,  and  an 
imaginary  part,  B. 

Therefore,  for  each  frequency  component. 


Subroutine  listings  are  given  for  REORDER  CFFTS,   CFFTRC,  and 
REALTRAN  on  pages  11-12,  pages  13-14,  pages  15-16,  and  page  17, 
respectively. 

These  four  subroutines  are  part  of  an  FFT  subroutine  package 
contained  in  the  public  user's  subroutine  library  of  the  in- house  time 
share  computing  system  at  NBS,  Boulder,  Colorado.     For  additional 
information  concerning  these  subroutines  see  Ref.  [5]. 


MAGNITUDE  (r) 


(2.  5) 


PHASE  (r) 


(2.6) 
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01  MENS  UN  A(1^24)>d(liJ24)*LSr  (lb)*LG(lb)*Sr(l:)) 
CJMMJN  M*LC*ST 
M  =  MM 

CALL  FFTC 

533  CJNriNJE 

IF(M)301 #173*331 
3131  CJNriNUE 

JB=M-1 

K8=3 

1=3 

KJ=LC( JA)-l 

0 J  23  KA  =   1 >KJ*2 

T   =  A<KA+l  ) 

A(KA+1 )   =  3 (KA) 

3(KA)   =  T 
23  CJNTINJE 

IF (M-1 )332> I 73>332 

332  CONTINUE 
LIM=M/2+l 

33         KS   =  LC( JB+1 )+Ka 
KJ  =  KS 
JJ=LC( JA-J8) 
KK  =  K3  +  JJ 
43         K  =  KK  +  JJ 
53  r   =  A  <KK  +  l  ) 

A(KK+1 )  =  A(KS+1 ) 
A(KS+1 )  =  T 
T  =  B(KK+1 ) 
B<KK  +  l )  =  B (KS  +  l  ) 
B<K3+1 )  =  T 
=  KK  \ 
KS  =  KS  +  1 
IF(KK-K)53*333*333 

333  CJNTINJE 
KK=KK+JJ 

KS   =  KS  +  JJ 
IF (KK-KU )40*334>334 

334  CJNTINJE 
IF(J8-LIM)93*93>335 
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306  CONTINUE 
J8=JB-1 
I   s  I-^l 
^  LST< I)   =  J3 

GO  TJ  3^ 
90  IF(I)17a*170#306 
306  CONTINUE 
JB=LST(I) 
I   -   I  -  I 
KB  =  KS 
GJ  TJ  30 
170  RETURN 
END 
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SJQ'RJJTIiME  CFFrS  (A*3*MM*SCALii:*iMc:XP) 
DIM£t\JSI  JiM  A(163  84)*a(16384)#  JC(l:>)>SNr<lb) 
CJMMJN  M>JC*SNT 
M=MM 
MA  =  MM 
CALL  FFTC 

IF(M-23i<3)53'3*l  73>5^3 
6  3  3  CJiMTINJE 
iM=JC(M+l  ) 
NH=iM/2 
NQ=NH/2 

SC   =  SCALE 
IF(M-2 )301 *332*3a2 

331  IF(M)123>123> 101 

332  CJNTIiNiJE 

EXPS   =   ISIGiMCl  *N£XP) 
NiM   =  N-1 
K   =  1 

03   100  JA  =  2*M 
C£=SNT<MA) 

I 

*  CE 


33 
43 


MA 

MA   -  I 

CD 

2.   *  CE 

SO 

-SNTCMA) 

R  = 

-2'   *  CD 

CN 

1  • 

CM 

0. 

S.M 

0. 

JJ 

0 

KK 

1 

SM 

■^EXPS 

JSPAN  =  NH 

JSPAiM/2 

KS 

KK  +  JSP 

RE 

s 

A(KK)  - 

A(KS) 

ACKK)   =  A(KK)  +  A(KS) 

FIM  «  8(KK)   -  a(KS) 

8CKK)   =  8(KK)  +  a<KS) 

A(KS)   =  CN*R£  -SN*FIM 

8<KS)   =  SN*R£  CN*FIM 
KK  =  KK  +  NH 
KS   =  KS   +  NH 
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R£   =  A(KK)    -  A(KS) 
A<KK)   s  A(KK)   +  A(KS) 
FIM  =  aCKK)  -8<KS> 
^  B<KK)   =  3(KK)   +  8(KS) 

8(KS)   =  SM  *  RE   +  CM  *  FIM 
KK  =  KS  NH 
IF<KK-N)40*333#3a3 
3  33  CONTINUE 
50         KK  s  KK  -  NN 
JJ  *  JJ  +  K 
IF( JJ-NQ)304*90#93 

304  CONTINUE 

60         CD  =  R  ♦  CN  +  CO 
CN  =  CD  CN 
SO   s  R  *  CM  +  SD 
CM  *  CM  +  SO 
SNa-CM*EXPS 
SM=CN*EXPS 
GO  TO  40 
90         K  s  K  +  K 
100  CONTINUE 
101         DO  1 10  KK  s  1 *N#2 
KS  =  KK  +  I 
RE  =  A(KK)    -  A(KS) 
ACKK)   s  (ACKK)   +  ACKS)) 

ACKS)   =  RE 
FIM  =  8CKK)   -  8CKS) 
8CKK)   =  C8CKK)   +  SCKS)) 

8CKS)   =  FIM 
110  CONTINUE 
120       IF  CABS  CSC -1 • )-l .E-1 3)13b*335>305 

305  CJNTINUE 

DO   130  Ja=l^N 
A<J3)   =  SC  *  ACJ3) 

BCJB)   =  SC         *  8CJ3) 
130  CONTINUE 
135  RETURN 
173  RETURN 
END 
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"FFTS" 

SJB^JtJTINE  GFrT^C    (  A  *  3  *  MM,  SCALd  *  ) 
131M£NSI  JN  A(l  /)^4)>3(l/)34)*  JlJCl  is  )  >  S  (  1  3  > 
GJMMJ.vJ  M*JD*S 
M  =  MM 

CALL  FFTC 

IF(M-2^^'4)bi3<3*17^>oa^ 

l^J  =  JL)(M■^l  ) 
K   =  iM/4 

=  K 
iMiM   =  iNj-l 
J3PAN  =  1 

SC   =  SCALE 
IF  CASS (SC- 1 . )-l  .£-1 ^)7*30l *3J1 
3(dl  CONTINUE 
6  OJ  i>  JC   =   1  *iM 

A(JC)   =  SC  *  A(JC) 

8<JC)  =  SC       *  a(jc) 

5  CJi^JTIi^aE 
7  IF (M)332* 170*302 

332  CJNriNUE 

OJ   13  KK  =   1  >iM,2 
KS   =  KK+1 

>^£   =  A(KK)    -  A(KS) 
A(KK)    =  A(KK)    +  A(KS) 

A(KS)   =  RE 
FIM        =  8(KK)    -  aCKS) 
8(KK)    =  B(r<K)   +  3(KS) 

aCKS)    =  FIM 
'   10  CJiMTIiMJE 

IF(M-1 )303*170,303 
303  CJNTIiMJE 

EXPS  =   ISIGiMCl  ,N£XP) 
i>iP  =  1 
DJ  90   J3   =  2*M 
SD   =  -S< Ja-1 ) 

CD  =  2.*  scja)  *  s(ja) 

R   =  -2»  *  CD 

CN  =  1  . 

CM  =  0. 

SlM  =  0» 

JJ    =  0 

KK   =  I 
SM  =  +EXPS 
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12  JSPAi>JH  =  JSPAN 

JSPAiM   =   JSPA.M   +  JS-^AiM 
23  KS   =  KK  +  JSPAiNJ 

RE  =  CiM  ♦  A<KS)    -  *  BCKS) 

FIM  =  Si^  *  A(KS)    +  Gi>J  *  B(KS) 
A(KS)    =  A(KK)    -  R£. 
A(KK)   =  A(KK)   +  RE 
3(KS)   =  aCKK)    -  FIM 
B(KK)   =  a(KK)  FIM 
KK   s  KK  +  JSPANH 
KS   =  KS   ♦  JSPANH 
FIM  =  SM  *  A(KS)   +  CM  *  B(KS) 
RE  *  CM  *  A(KS)    -  SM  *  t3<KS) 
ACKS)   =  A(KK)    -  RE 
ACKK)   =  ACKK)   +  RE 
3CKS)    =  3CKK)    -  FIM 
3CKK)   =  3CKK)    +  FIM 
KK  «  KS   +  JSPAiNJH 
IFCKK-N)2a*304*304 
334  CONTINUE 

30  KK  =  KK  -  NN 

JJ  =  JJ  ■•■  K 
IFC JJ-NQ)305*80>80 
305  CJNTINUE 
40         CD   s       *  CN  +  CD 
CN  =  CD  +  CN 
SM=CN*EXPS 
SD  =  R   ♦  CM  +  SD 
CM  =  SD  +  CM 
SN=-CM*£XPS 

GJ  TO  20 
30  K  s  K/2 

RAO=.b*RAD 
90  CONTINUE 
170  RETURN 
END 
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SUBRJLiriNE  REALTRAN(A#9>MM*;M£*  IN^) 

DIMEiMSIJi^  ACia24)*B(ia2  4)*JC<lb)*SrClb) 
CJMMJiM  M*JC*ST 

1  MsMM 
CALL  FFTC 

IF(M-20^  )i>30>l  70*b3i3 
500  CJNTINUe 

K=JC(M+l) 

N«K 
iMH=N/2 
NK=MH+l 

CN  =  ISIGf^(l  >  IN\/) 
S»^=ISIGiN)(l  *iM£) 
IF(CN)301 #3^3 

301  CONTINUE 

2  f IM=A(1  )-A<N+l  ) 
Ad  )=A(1  )+A(N+l  ) 
8(1 )=FIM 

GJ  TJ  4 

3  A<N+1 )s2.*(A(l )-B(l ) ) 
ACl  )«2.*(A<1 )+B(l ) ) 
B(N+1 )=0 

8(1  )=0 

4  IF(M)302*170*302 

302  CJNTINUE 
SD*CN*SN*ST (M) 
F?a2  .*sr  (M+l  ) 

R=-R*R 

CD=-.5*CN*R 

SN=0. 

DJ  5  J=:2*NK 

CD   =  R   *  CN  +  CD 

CN  =  CD  +  CN 

SO   =  R   *  SN  +  SD 

SN  =  SD  +  SN 

AA  =  A (J)   +  A(K) 

A8   =  A(J)    -  A(K) 

BA=3(J)+8(K) 

3B   =  B(J)    -  B(K) 

RE   =  CN  *  BA  +  SN  ♦  AB 

FIM  =  SN  *  BA  -  CN  *  AB 

B(K)   =  FIM  -  BB 

8(J)   =  FIM  +  88 

A(K)   =  AA   -  RE 

A(J)   =  AA   +  RE 

5  K=K-1 

170  RETURN 
END 


17 


Two  FORTRAN  II  programs  are  described  in  this  section.  They 
are  used  to  modify  the  time  domain  waveform  data  stored  in  local 
memory,  making  these  data  suitable  for  frequency  domain  transfor- 
mation. 

The  first  program,   stored  in  a  user's  file  in  the  time  share 
system  under  file  name  /XXI/,  performs  three  functions.  However, 
before  a  detailed  description  of  these  functions  is  presented,  it  is 
necessa^ry  to  discuss  the  format  of  the  data  retrieval  from  local 
memory. 

The  time  domain  waveform  data  is  stored  locally  in  a  1000  point 
BCD  memory.    Each  Y-axis  data  point  consists  of  three  decimal  digits. 
The  X-axis  data  is  not  needed  since  the  Y  data  is  stored  and  trans- 
mitted sequentially.    This  1000  point  waveform  is  converted  to  ASCII 
teletype  code  and  put  on  paper  tape.    It  is  then  read  into  the  time  share 
system  from  a  paper  tape  reader  and  placed  in  file  name  /RAW/, 
Program  /XXI/  is  then  implemented  as  follows: 

1.  The  first  f\inction  of  /XXI/  is  to  shift  the  array  indexing 
of  /RAW/,    An  idiosyncracy  of  the  hardware  is  that  the 
first  Y  data  point  appears  on  paper  tape  three  times  in 
succession.     Thus  it  is  necessary  to  delete  the  first  two 
data  points  in  file  /RAW/.    This  is  accomplished  by  index 
shifting  such  that  M(1)  =  M(3),  M(2)  =  M(4),  etc, 

2,  At  this  point,  /RAW/  consists  of  1000  three-digit  Y  data 

points  representing  a  time  domain  waveform.     The  FFT 

algorithm  requires  the  waveform  data  array  to  be  of  length 
N 

2     ,  where   N    is  an  integer.     The  closest  length  is  1024 
data  points  so  24  points  must  be  added  to  the  array.  This 
may  be  accomplished  in  any  suitable  manner,  but  in  order 
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to  introduce  a  minimum  of  error  to  the  data,  the  following 
method  has  been  chosen.    The  average  value  of  the  last  ten 
data  points  is  computed.    This  value  is  then  used  to  fill  the 
array  by  placing  it  in  positions  1001  thru  1024.    In  other 
words,  the  averaged  final  value  of  the  waveform  is  simply 
repeated  24  times. 
3.      The  third  program  function  is  a  form  of  signal  averaging. 
This  portion  of  the  program  is  optional.    If  an  unaveraged 
waveform  array  is  desired  it  is  stored  in  file  name  /DON/. 
If  the  averaging  process  is  required,  the  data  output  appears 
in  file  name  /DONER/.     Waveform  averaging  is  accomplished 
by  comparing  the  value  of  each  data  point  with  the  average 
value  of  its  neighbors  on  either  side.    If  the  value  of  the 
tested  data  point  exceeds  certain  limits,  that  value  is  changed 
to  equal  the  neighbors'  averaged  value.    See  Figure  2.  1. 

Program  /XXI/  receives  values  for  NR  and  NL  from  the  TTY 
keyboard  at  run  time.    From  Figure  2.  1  it  is  apparent  that  NR  is  the 
total  number  of  neighbors  nearest  the  test  point.    In  this  example  NR  = 
12,  or  six  on  either  side  of  the  test  point.    LOAV  is  the  computed  aver- 
age value  of  the  six  data  points  to  the  left  of  X    ,  and  LIAV,  the  average 
of  the  six  points  to  the  right  of  X  .     The  average  value  of  LOAV  and  LIAV 
then  is  computed  and  stored  as  LOT.     Therefore,  LOT  is  the  average 
value  of  the  NR  nearest  neighbors  of  X  .    The  absolute  value  of  X  -LOT 
is  then  compared  with  the  value  chosen  for  NL.    If    Ix^-LOT  |  exceeds 
NL,  then  X^  is  reset  equal  to  LOT.    If    |X^-LOT|  is  less  than  or  equal 
to  NL,  then  X    is  unchanged.    NR  must  be  a  positive  even  integer  and 
NL,  a  positive  integer. 

Figure  2-2  shows  a  typical  time  domain  waveform  along  with  the 
results  of  averaging  for  various  values  of  NR  and  NL  using  Program 
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/XXI/     The  program  listing  is  as  follows: 

NOTE:    Because  a  detailed  explanation  of  both  the  input  data  and  the  out- 
put data  programs  is  contained  in  the  text,   comment  statements 
are  not  included  in  the  program  listings. 
/XXI/ 

r  DIM£i>4SIJ>J  M(1024) 

QIMEi^JSIJN  L(1324) 
JPEN   (3*  INPJT*/i^AW/) 
JPEN  (4*0jrPJr*/0JNER/) 
READ  3* 1 3* (M(K) *K=l *  1 3 33 ) 
TYPE  25 

25  FJRMAT    (10HNL  AND  NR ? ) 

ACCEPT  lb*NL*NR 
15  FORMAT  (212) 

0  J  5    1  =  1 *1000 
:>  i-U  I  )=i^(  1+2  ) 

D  J  7    1  =  991  »\d'dd 

7  SJM=SJM+M(I) 
A\/G=SJM/10 

OJ   8  1=1^31*1324 

8  M(I)=Ay/G 
CLJSE  (3) 

I  JPEN   (5*  Jjr?>JT*/D  JiM/) 

DJ   17  J=l*1324 
17  WRITE  5*23*(M(J)) 

13  FJRMAT  (1333(13/)) 

23  FJRMAT  (1324(13/)) 

0 J  35  K  =  l *  1324 
35  L(K)=M(K) 

DJ  4  3   J=l * ( 1 324-NR ) 

LJ=3 

LI  =3 

OJ  60  K  =  J*  (  J+(.MR/2  )-l  ) 
60  LJ=LJ+L(K) 

OJ   73   K=( J+(NR/2 )+l ) * J+NR 
73  LI=LI+L(K) 

LJA\/=LJ/  (NR/2  ) 

LIA\/  =  LI/(i>JR/2  ) 

LJT  =  (LJA\/+LIA\/)/2 

IF ( IA8S (L( J+(NR/2 ) ) -L JT ) -NL) 43  *  43  * 2 1 
21  L(  J+(iMR/2  )  )  =  LJT 

43  CJiMTINJE 

OJ  83  J=W1324 
83  WRITE  4*23* (L( J) ) 

ST  JP 

END 
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The  second  data  input  program  is  called  /MODI/.     The  input 
data  to  this  program  is  listed  in  file  name  /DONER/  and  is  the  output 
from  program  /XXI/.    Since  the  FFT  algorithm  will  produce  erroneous 
results  for  non-periodic  waveform  inputs,  this  program  insures  that  the 
input  waveform  will  appear  periodic.     This  is  accomplished  by  first 
deleting  all  the  odd-indexed  data  points  in  the  input  waveform.  Then 
the  resulting  512  point  waveform  is  inverted,  vertically  shifted,  and 
placed  in  indices  513  through  1024.    See  Figure  2.  3  for  a  sample  wave- 
form both  before  and  after  being  made  periodic. 

The  output  of  this  program  is  then  in  proper  form  to  be  trans- 
formed using  the  FFT  subroutines.  The  /MODI/  listing  is  as  follows: 


/MJOl/ 

DIMENSION  K<1024) 
DIMENSION  L(1024) 
DIMENSION  M<ia24) 
DIMENSION  N<1024) 
OPEN  <3*  INPUT  */OONE!^/) 
READ  3*10* CMC J)# Jsl > I  324) 

10  FORMAT  CI3) 

DO  5  J=l*5l2 

5  NCJ)=:M<J*2) 
ISJM=0 

DO  20  J=l *  10 
20  IS'JM=IS'JM-»-N(  J) 

IA^G=ISUM/10 
•  ISHIFT  =  N(5l  I  )-IA\/G 

DO  30  J=l *512 
30  L(J)=N(J) 

DO  40  J=l *512 

40  K< J)=2*N(51 1 )-L( J)-lSHIFr 

DO  50  J=l *512 
50  N( J+512)=KC J) 

OPZi^   (4*OJTPUT*/O0N£/ ) 

00  70  J=l *  1 024 
70  WRITE  4*60*N<J) 

60  FORMAT  (13/) 

STOP 

END 
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2.  3    Output  Data  Programs 

The  first  data  output  prograna  is  called  /TRI/.     This  program 
accepts  the  output  of  /MODI/,  a  1024  point  averaged,  periodic  time- 
domain  waveform  stored  in  file  name  /DONE/  and  performs  the  FFT  of 
this  waveform.     The  data  output  from  this  program  is  available  in  the 
following  files: 

/DONEA/  =  A(r) 
/DONEB/    =  B(r) 

/DONEX/    =  I  [A(r)]^  +[B(r)]^  .^^^ 

where  S(r)   =  A(r)  +  jB(r) 

and  S(r)  are  the  512  xinscaled  complex  frequency  coefficients  associated 
with  the  1024  point  input  waveform.    /DONEX/  is  really  a  scaled  dis- 
creet magnitude  spectrxxm  of  the  input  waveform.    Scaling  is  accomplished 
by  multiplying  each  A(r)  and  B(r)  (with  the  exception  of  A(0))  by  2/N  where 
N  is  the  length  of  the  transform.    In  this  case  N  =  1024.    A(0),  or  the 
D.  C.  term,  must  be  multiplied  by  1/N  for  proper  scaling. 

The  second  data  output  program  is  called  /LSP/  and  is  used  to 
scale  the  output  magnitude  spectrum  /DONEX/  so  it  may  be  returned  to 
local  memory  and  displayed  on  a  CRT.     The  scaling  equation  used  to  yield 
an  integer  spectrum  in  the  range      0  ^  M(r)  <  1000  is: 

M(r)   =  B'(r)  =   200  1ogjQ[l00  |S(r)|] 

where    |S(r)|    are  the  scaled  magnitudes  of  the  complex  frequency  co- 
efficients, B  (r)  are  the  rescaled  decibel  values  of    |S(r)|,  andM(r)  are 
the  nearest  integer  values  of  B'(r).    See  Figure  2.  4  for  some  sample 
displays  of  magnitude  spectra  in  dB  versus  frequency. 
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The  third  data  output  program,  called  /UNTRI/,  is  used  to  per- 
form the  inverse  FFT  on  frequency  domain  data.    The  inputs  to  this 
program  are  the  unsealed  complex  frequency  coefficients  from  /TRI/. 
They  are  contained  in  file  names  /DONEA/  and  /DONEB/.     The  outputs 
of  this  program,  contained  in  /Yl/  and  /Y2/,  are  the  first  and  last  half, 
respectively,  of  the  time  domain  waveform  associated  with  /DONEA/  and 
/DONEB/.    The  /UNTRI/  program  was  used  only  to  check  the  validity  of 
the  FFT  program  package.    The  program  listings  for  /TRI/,  /UNTRI/ 
and  /LSP/  are  as  follows: 


OIMiNSIJN  MC1324)>A(513),"3(513)>C(1.3:'4).0(513) 

R£AD  3>13* CM<K)^K=1 > 1024) 
13  FORMAT  (13) 

OJ  43  J=l > I  324 
40  C(J)=M(J) 

DJ  53  J=l ,512 
53  A<J)=C<J) 

DO  60  J=l  *512 
63  8CJ)=C< J+512) 

30  FORMAT   <F1 1  .2) 

CALL  REORDER  (A>a#9) 

CALL  CFFTRC   <A,3*9# .5*1 .* 1 ) 

CALL  REALTRAN  CA*8^9*1,1) 

CLOSE  (3) 

OPEN  <4>0UrPUT> /DONEA/) 
WRITE  4*33* (A(J)> J=l >513) 
OPEN  <5*0UrPUT>/00NEa/) 
WRITE  5>33* (3(J), J=l *513) 
DO   100  J=l,513 
133         DC J)  =  CSQRT<  CA<J)**2)  +  (8<J)**P) ))/5iP, 
CLOSE  (4) 
CLOSE  <5) 

OPEN   (6,0UTPUT^/00NEX/)  • 
WRITE  6*120> (DCJ)* J=l >513) 
123         FJRMAT  (F7.2) 
STOP 
END 
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OlMENSIJiNJ  A<5l3):.3(bl3) 
J-^EN   C3  >  INPJr^/D  JiMEA/) 
J-^EN   (4>  lNPJT*/0  JiNJES/) 

^EAO  4*1 3> (8CK)*K=1 *513) 
13  FJ-Ri-JAr  CF11.3) 

CALL  REALr.^AiM   (A*B*  9> -1  , -I  ) 
CALL  CFFrS    CA>8# 9^ 1 ./ 1 354. > -I ) 

CALL  REJ^OEr^  (A*3*9) 
CLOSE  (3) 
CLJSE  (4) 

J-'EiM   (5> JUrPUT^/yi /) 
JPEN   (6* JjrPUr*/Y3/) 
W^ITE  5>1 0> (A<K)>K=1 *D13 ) 
w^irE  6>l'a*  (a(K)jK=l^:5l3) 
STJP 
END 


NOTE:    In  the  listing  which  follows  B(K)  represents  the  term 
B'(r)  as  defined  on  page  22. 


/LS-^/ 

0  I  MEiNJS  1        A  (  5  1  3  )  *  3  (  1      4  )  *  H  (  1  J  a  4  ) 

OPEN   <3*  IiMPiJT>/D  JiMEX/ ) 

READ  3*l'^>  CA<K)*K=1  *:>13) 

A(l)=3. 

OJ   15  K=l >256 
15  AC2*K)=2k)0.*<ALJG  (  1  3<3  .*A  (2*K)  )  )/ CALJGC  1  J  .  J)  ) 

D J  25   K=l ^512 

3 (2*K-l )=0. 
25  8C2*K)=ACK) 

DO  35  K=l >1024 
35  M(K)=8CK) 

OPEN  <4*0UTPUT>/SPEC/) 

WRITE  4#22l*  <M(K)*K=sl  ^  1  324) 

OPE,i^  (5*0UTPLfT*/SHJRr/) 

WRITE  5*20* CM<4*K)*K=1 *256) 
10  FORMAT  <F7.a) 

20  FORMAT  <I3) 

STOP 

END  '  ■  ■ 
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3.    EXPERIMENTAL  CONCLUSIONS 

3.  1    Problems  Encountered  Using  the  FFT 

The  primary  goal  of  the  work  upon  which  this  report  is  based  is  to 
develop  both  hardware  and  software  for  a  particular  system.    At  this 
point  in  time,  the  main  desired  function  of  the  system  is  to  produce  the 
gain  versus  frequency  characteristics  of  certain  electrical  networks 
using  time  domain  techniques.    As  a  result,  only  specific  FFT  applica- 
tion problems,  pertinent  to  this  particular  system,  will  be  discussed. 
Problems  of  a  more  theoretical  nature  are  discussed  in  detail  in  the 
literature,  e.g.,  [3]  or  [4]. 

Basically,  the  system  works  in  the  following  manner.    A  small 
transition- time  voltage  step-waveform  is  generated,  observed  with  a 
sampling  oscilloscope,  and  the  oscilloscope  analog  output  recorded. 
The  FFT  of  the  recorded  waveform  is  then  computed,  and  a  magnitude 
versus  frequency  graph  is  recorded.     The  voltage  step-waveform  is  then 
applied  to  the  input  terminals  of  the  network  under  test  (e.  g.  ,  a  coaxial 
attenuator);  the  network  output  waveform  is  observed  with  the  sampling 
oscilloscope,   recorded  and  transformed.    The  logarithmic  difference 
of  these  two  transforms,  then,  is  computed  to  yield  the  gain  versus 
frequency  characteristics  of  the  network. 

Actually,  four  major  problems  or  limitations  were  encountered 
which  had  a  detrimental  effect  on  network  measurements.    Each  problem 
shall  be  discussed  individually. 

The  first  problem  was  due  to  the  choice  of  a  voltage  step  as  an 
input  waveform.     The  FFT  algorithm  assximes  the  input  waveform  to  be 
periodic.    Thus  an  unwanted  discontinuity  will  appear  between  the 
beginning  and  end  points  of  the  step  waveform,  implying  that  the  voltage 
step  was  instantaneously  turned  off.    The  FFT  will  transform  the  in- 
stantaneous transition  as  well  as  the  valid  voltage  step  waveform  data, 
yielding  a  grossly  invalid  frequency  spectrum. 
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The  objective,  then,  was  to  make  the  step  waveform  appear 
periodic  while  introducing  neither  discontinuities  nor  other  non- causal 
voltage  changes.     This  was  accomplished  through  the  use  of  FORTRAN 
-  II  program  /MOD  1/.    This  program  modifies  the  input  step  waveform 
in  such  a  manner  that  it  appears  to  the  FFT  as  if  the  voltage  step  were 
turned  off  in  exactly  the  same  manner  in  which  it  was  turned  on.  (See 
Figure  2.3).    Consequently,  the  modified  waveform  appears  periodic 
with  no  discontinuities  and  no  new  non-causal  data  added  to  it. 

The  second  problem  was  also  solved  using  program  /MOD  1/. 
Figure  3.  1  consists  of  two  photographs  of  the  magnitude  spectra  in  dB 
of  two  ideal  rectangular  pulses  of  arbitrary  width.     These  spectra  are 
both  difficult  to  interpret  since  they  are  made  up  of  two  curves  super- 
imposed upon  one  another.    In  fact,  one  curve  on  each  photograph 
represents  the  magnitude  spectr\im  of  the  even  harmonics  while  the 
other  curve  represents  the  magnitude  spectrum  of  the  odd  harmonics. 

It  was  found  that  the  manner  in  which  program  /MOD  1/  produces 
a  periodic  waveform  from  a  voltage  step  waveform  also  ensures  that 
the  even  harmonic  magnitude  spectra  will  be  forced  to  zero  and  the 
odd  harmonic  magnitude  spectra  will  decrease  monotonically  as  in 
Figure  3.  2.    In  terms  of  Fourier  series  analysis,  program  /MOD  1/ 
generates  a  single  period  of  a  periodic  waveform  in  which 

f(t)   =    -  f  (t  +  ^  )  (3.  1) 

where    T   is  the  period  of  the  program- generated  waveform.  From 
Fourier  series  analysis,  any  periodic  waveform  that  fulfills  equation 
(3.  1)  will  contain  only  odd  harmonics. 

The  third  problem  is  really  a  limitation  inherent  in  the  use  of  a 
step  waveform  as  a  driving  function.     The  step  waveform  must  be  made 
periodic  in  order  to  satisfy  the  requirements  of  the  FFT.     Thus,  the 
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modified  waveform  will  always  resemble  a  rectangular  pulse.    But  the 
Fourier  integral  transform  of  a  rectangular  pulse  is  of  the  form 
six  x/x  where   x   is  a  function  of  both  frequency  and  the  width  of  the 
time  window.    Consequently,  there  is  an   inherent    1/x   decay  of 
signal  input  level  as  a  function  of  frequency.    This  signal  level  decay 
may  be  overcome  somewhat  by  careful  selection  of  the  time  window 
used,  and  possibly  by  pulsed  RF  techniques. 

The  fourth  major  problem  is  that  of  noise.     This  problem  is  not 
caused  by  use  of  the  FFT,  but  rather,  it  is  caused  by  the  data  acquisition 
system.    It  is  worth  mentioning    because  several  software  problems 
result  in  processing  the  resultant  noisy  signals.     The  magnitude  spectra 
shown  in  Figure  2.  4  are  noisy  spectra  typically  encountered.    As  may 
be  seen,  only  the  very  beginning  of  each  spectrum  appears  to  show  any 
continuity.     The  rest  is  much  too  noisy  to  be  interpreted. 

Aside  from  attempts  to  diminish  the  noise  sources  within  the  hard- 
ware, two  software  techniques  were  used  to  bring  the  noise  level  down. 
One  was  signal  averaging  as  in  FORTRAN  II  program  /XXI/.     The  other 

was  taking  transforms  over  several  decreasingly  smaller  time  windows 
to  obtain  usable  frequency  domain  data  over  an  increasingly  wider 
frequency  band. 

3.  2    Time  Domain  Measurement  Example 

In  this  section  the  measurement  of  the  gain  versus  frequency 
characteristics  of  a  coaxial  10  MHz  low-pass  filter  is  presented.  Figure 
3.  3  is  a  block  diagram  of  the  waveform  data  acquisition  system  used  to 
perform  this  measurement.    For  purposes  of  this  report,  this  data 
system  simply  acquires  and  stores  in  local  memory  a  1000  point  repre- 
sentation of  a  time  domain  waveform.    The  input  step  waveform  and  the 
low-pass  filter  output  waveform  are  shown  in  Figure  3.  4.    In  the  example, 
the  input  waveform  is  unsmoothed  while  the  output  waveform  has  been 
smoothed  using  FORTRAN  II  program  /XXI/. 

27 


The  FFT  software  programs  were  then  applied  to  both  waveforms 
and  the  logarithmic  magnitude  versus  frequency  spectrum  of  each  wave- 
form was  recorded  in  local  memory.    An  X-Y  recording  of  these  two 
spectra  is  shown  in  Figure  3.  5.    This  figure  clearly  shows  the  1/x 
decay  characteristics  of  the  transformed  step  waveform.    Also,  it 
exhibits  a  relatively  noiseless  gain  window  of  about  60  dB  below  the 
fvmdamental  frequency  amplitude. 

Figure  3.  6  is  an  X-Y  recording  of  the  algebraic  difference  of  the 
input  and  output  spectra  shown  in  Figure  3.  5.     This  figure  shows  a  fairly 
flat  response  of  the  low -pass  filter  out  to  approximately  9  MHz.  The 
gain  then  drops  at  a  rate  of  about  20  dB/MHz  into  the  noise. 

The  results  of  this  measurement  appear  to  be  promising.  The 
purpose  of  this  measurement  was  merely  a  qualitative  check  of  the  basic 
system  hardware  and  software,  but  the  resulting  curve  of  Figure  3.  6 
shows  a  noise  level  and  a  frequency  resolution  that  indicates  the  feasi- 
bility of  the  time- domain  calibration  system. 

3.  3    Future  Developmental  Considerations 

With  the  goal  in  mind  of  developing  a  time  domain  network  calibra- 
tion system  competitive  with  current  capabilities  of  frequency  domain 
calibration  systems,  much  more  work  is  necessary.    Outlined  below  is 
a  schedule  of  tasks  deemed  necessary  to  meet  this  goal. 

1.  The  first  task,  that  of  implementing  a  dedicated  minicomputer 
into  the  system,  is  already  underway.    The  minicomputer  will 
allow  much  higher  speeds  of  data  handling  and  processing  as 
well  as  increased  system  automation  capabilities. 

2.  After  implementation  of  the  minicomputer,  a  more  detailed 
and  exhaustive  set  of  experiments  aimed  at  error  analysis 
may  be  performed.    At  present,  use  of  the  time  share  com- 
puting system  is  much  too  time  consuming  to  permit  an  in- 
depth  study  of  systematic  or  random  errors. 

28 


3,  With  sources  of  errors  pinpointed,   software  may  then  be 
developed  as  far  as  possible  both  to  determine  the  magnitude 
of  and  to  reduce  the  effects  of  those  errors. 

4.  Concurrent  with  software  development,  the  fourth  task  would 
be  hardware  development.     This  would  involve  development 
of  fast- rise  noiseless  step  generators,  improved  sampling 
systems,  increased  data  acquisition  resolution,  and  increased 
system  automation. 

The  general  concept  of  using  the  FFT  to  compute  the  gain  versus 
frequency  characteristics  of  electrical  networks  appears  quite  feasible. 
Most  probably,  the  real  limitations  on  measurement  precision  and 
accuracy  lie  in  system  hardware,  and  not  in  FFT  programming  or  other 
software. 

This  report  has  been  concerned  with  three  main  topics:    a  definition 
of  the  fast  Fourier  transform;  a  description  of  software  programs  used  to 
compute  the  FFT  on  the  in-house  time  share  computing  system;  and  a 
brief  description  of  problems  encountered  in  the  data  reduction  process 
used. 

The  next  report  on  this  continuing  work  will  be  addressed  to  the 
actual  use  of  the  system  for  quantitative  measurements  of  amplitude  and 
phase  characteristics  for  microwave  attenuators,  filters,  and  mismatches. 
These  measurements  will  also  be  compared  to  those  made  at  other 
laboratories  using  time  domain  techniques  (e,  g,  ,   see  Ref.  [6]  and  [7]  ). 


29 


REFERENCES 

1.  Cooley,  J.  W.  and  Tvikey,  J.  W.  ,  "An  Algorithm  for  the  Machine 
Calculation  of  Complex  Fourier  Series",  Mathematics  of  Com- 
puters, Vol.  19.,  pp.  297-301,  April,  1965. 

2.  Cochran,  W.  T.  ,  et  al.  ,  "What  is  the  Fast  Fourier  Transform?  ", 
IEEE  Trans,  on  Audio  and  Electroacoustics,  Vol.  AU- 15,  No.  2, 
June,  1967. 

3.  Cooley,  J.  W.  ,  et  al. ,  "Application  of  the  Fast  Fourier  Trans- 
form to  Computation  of  Fourier  Integrals,  Fourier  Series,  and 
Convolution  Integrals",  IEEE  Trans,  on  Audio  and  Electro- 
acoustics,  Vol.  AU-15.  No.  2,  pp.  79-84,  June  1967. 

4.  Bergland,  G.  D.  ,  "A  Guided  Tour  of  the  Fast  Fourier  Trans- 
form", IEEE  Spectrum,  pp.  41-51,  July,  1969. 

5.  FFT  subroutine  package  written  by  David  L.  Lewis,  National 
Oceanic  and  Atmospheric  Administration,  Space  Environment 
Laboratory,  Boulder,  Colorado  80302. 

6.  Nicholson,  M.  A.  ,  et  al.  ,  "Applications  of  Time-Domain  Metrology 
to  the  Automation  of  Broad-Band  Microwave  Measurements",  IEEE 
Trans,  on  Microwave  Theory  and  Techniques,  Vol.  MTT-20,  No.  1, 
January  197  2. 

7.  Cronson,  H.  M.  ,  et  al.  ,  "Time  Domain  Measurement  Study", 
Sperry  Rand  Research  Center,  Sudbury,  Massachusetts,  August, 
1972. 


30 


in 

CO 


o" 

II 


•■-I 

u 
> 


u 
o 

(U 

> 

M-t 

o 

cn 
•+-> 
I— I 
:i 
t» 
<u 

CO 

O 

> 
a 


u 
o 

> 


o 


CO 
.1-1 
O 

1 — I 

a 
u 
•i-i 


bio 


32 


34 


cr 


o 

a 
u 

o 

<u 

to 
o 

?J 
<u 
?s 
cr 
<i) 

■<+^ 

CO 

> 

<i; 
lis 

-I-l 

fl 

two 


o 

01 
CO 

•  r-l 
t— I 

d 
o 


•CO 


*0 
-1-1 

F=4 


M 

O 
m 

> 
.1-1 

a 

o 


35 


37 


W 

H 

OMPU 

J 


o 

I— I  >- 
I—  UJ  q; 
llj  ct:  o 

r—   ^  O  2: 

C3  O  LjJ 


LlJ 


\  U-l 


O 
O 


Q  LU 


CD  Q 
1—1  o 
Q  CQ 


0 

0 

00 

0 

1 — 1 

 1 

UJ 

cc: 

LjJ 

0 

ai 

1— 

^  1— 

cc: 

1— f  1 — ( 

,  Q 

(— 

Q-  ZD 

Q 

<c  0 

<: 

CO 

00  1— 1 

X 

0 

_i  Q  q: 

CQ  UJ  UJ 

<:  >-  o 

•— '  ct  o 

q;  _i  >— I 

ct  UJ  q: 

>  Q  I— 


C3 
I — I 

_l 

CO 


UJ 
X  UJ 

00 


CO 


PE 

CD 

0 

0 

1 — 1 

LT) 

_J 

Q_ 

_l 

2: 

_l 

<c 

1 — 1 

00 

OSC 

Q 
O 


CD 
CD 


UJ  >- 

UJ 

Q  _J 

Q 

0  a. 

0 

1— 1  Q_ 

1 — 1 

Q  ID 

Q 

00 

_l 

_J 

UJ  00 

UJ 

2: 

Z  1— 1 

ZD  CQ 

ZD 

1— 

1— 

0  on 

u_  0 

UJ  h- 

3 

q; 

>  ■< 

0 

0 

2: 

3 

3  UJ 

1— 

UJ 

•  UJ 

U.  CD 
UJ 

CD 

cc: 

_J 

( — 1 

UJ 

<c 

UJ 

_i 

CD 

CJ 

_i 

Q- 

CD 

t— t 

Q- 

s: 

1 — 1 

1— 

2: 

<: 

Ql 

< 

00 

1— 

UJ 

1 

f 

> 

CD 

1 — 1 

Q 

_J 

<C 

Q_ 

UJ 

2: 

n: 

<c 

00 

38 


AtXI  00^ 


39 


41 


APPENDIX  A 

To  yield  some  insight  into  the  mechanisms  of  the  FFT,  an  eight- 
point  DFT  will  be  solved  below  using  the  FFT  algorithm. 

Given  a  time  domain  waveform  represented  by  eight  equally  time- 
spaced  data  points,  as  illustrated  in  Figure  A-1,  we  wish  to  compute 
the  DFT  for  this  waveform.    In  particular,  we  wish  to  calculate  the 

eight  A  's  associated  with  the  eight  assumed  X,  's.    The  DFT  for  the 
r  k 

set  of  data  points  is. 


:=o 


2  TT  jrk 

8 


r   =   0,  1,  2,  .  .  .  ,  7 


(A-1) 


For  the  first  reduction  of  transform  length  we  let 


B     =  >     Y,  e 

r  k 


8 


r   =   0,  1,  2,  3 


{A-2) 


and 


k=o 


■  4  ff  j  rk 

8 


r   =   0,  1,  2,  3 


(A-3) 


where 


and 


k  2k 


\  -  ^2k+l 


(A-4) 
(A-5) 


Using  the  simplifying  notation  that 

-2jrl 

N 


W    =  e 


we  have 


(A-6) 


.  =  >  X,  W 
r  k 


rk 


r   =   0,   1,  2,  .  .  . ,  7 


(A-7) 
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^  2rk 


and  B      =^  Y^W  ^  r   =    0,   1,   2,  3.  (A-8) 

k=o 


=J2  Zj^^^""^  r    =    0,    1,   2,   3.  (A-9) 

k=o 

Thus,  we  have  now  reduced  the  problem  to  that  of  finding  two  four-point 

DFT's  rather  than  one  eight-point  DFT.    (See  Figure  A- 2). 

For  the  second  reduction,    B      and   C      may  each  be  reduced  to 

r  r 

two  two-point  transforms  as  follows: 

Let  ^ 

D  R  W^^^  r    =    0,   1  (A-10) 

k=o 


and  E     =  VsW^^^  r   =    0,   1  (A- 11) 

r  k 
k=o 

where  R^^  =  Y^^  =  X^^  (A-U) 
and  \   =   ^2k+l    =  ^4k+2  '^"l^' 


and  let 


4rk 

F     =  V       W  r   =    0,   1  (A-14) 

r  ^ 


"Z)  ^^^""^  r    =    0,   1  (A-15) 

k=o 

where        T^^  =   Z^^  =  X^^^^  (A-16) 
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Then  D     is  the  two  point  transform  of  data  points   X     and   X^,    E  , 
r  o  4  r 

the  transform  of  X^  and  X.,     F  ,  that  of  X,  and  X^,  and  G  ,  that  of 

2  6        r  15  r 

X^  and  X^. 
3  7 

The  third  reduction,  which  would  yield  eight  one-point  trans- 
forms, is  really  unnecessary  for  the  DFT  of  a  single  point  function  is 
the  sample  itself.     Therefore,  substitution  of  the  data  points  into  (A- 10), 
(A- 11),  (A- 14),  and  (A- 15)  yields  the  following  set  of  equations. 


D  =    X    +W°X  F      =    X,  +W°X 

0  o  4  o  1  5 

D,  =    X    -  W°X,  F,    =   X    -  W°X^  (A-18) 

1  o  4  115 


E      =   X^  +W°X,  G     =   X^  +W°X^ 

o  2  6  o  3  7 


As  was  shown  in  the  text, 


and 


A     =   B    +W^C  r   =   0,   1,  2,  3  (A-19) 

r  r  r 


A,     ^.  =   B    -  W^C  r    =    0,   1,  2,  3  (A-20) 


r  r 


Likewise,  it  may  be  shown  that 


B     =  D    +W^^E  r   =   0,  1  (A-21) 


r  r  r 
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and 


=   D    -  W^'^E 
r  r 


C      =    F  +W^^G 
r  r  r 


Recalling  that 


W     =     -  W 


h) 


r  =  0,  1 


r    =    0,  1 


r        0,  1 


{A-22) 

(A-23) 
{A-24) 


(A-25) 


we  may  substitute  (A- 18)  into  (A-21)  and  (A-22)  and  then  into  (A- 19) 
and  (A- 20)  to  yield  the  following  matrix  equation: 


A 

o 

1  1 

1 

1 

1 

1 

1 

1 

X 

o 

^1 

1 

. 

1 

^1 

1 

-  1 

1 

-  1 

^3 

1 

+ 

-  1 

-w' 

^3 

A 

4 

1  -1 

+  1 

-  1 

+  1 

-  1 

+  1 

- 1 

X 

4 

^5 

1  -w' 

-  1 

^6 

1  -W^ 

-1 

1 

-1 

^7 

1  -w^ 

-  1 

+ 

3 



(A-26) 
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Examining  (A- 26)  we  see  that  half  of  the  multiplying  matrix  terms  are 
vinity.    In  addition,  the  matrix  is  symmetrical  about  the  main  diagonal. 
Consequently,  about  three -fourths  of  the  multiplications  are  unneces- 
sary.    Figure  A-3  is  a  signal  flow  graph  illustrating  this  computation. 
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Fig,  A- 2      Two  four-point  representations  of  the  same  time  domain 
waveform. 
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